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ABSTRACT 

Simulated dark matter profiles are often modelled as a 'NFW density profile 
rather than a single power law. Recently, attention has turned to the rather rig- 
orous power-law behaviour exhibited by the 'pseudo phase-space density' of the 
dark matter halo, which is defined dimensionally in terms of the local density and 
velocity dispersion of the dark matter particles. The non-power-law behaviour of 
the density profile is generally taken to exclude simple scale-free, in-fall models; 
however the power-law behaviour of the 'pseudo-density' is a counter indication. 
We argue in this paper that both behaviours may be at least qualitatively under- 
stood in terms of a dynamically evolving self-similarity, rather than the form for 
self-similar infall that is fixed by cosmological initial conditions. The evolution is 
likely due to collective relaxation such as that provided by the radial-orbit insta- 
bility on large scales. We deduce, from a distribution function given by first order 
coarse-graining, both the NFW-type density profile and the power-law pseudo- 
density profile. The results are not greatly sensitive to variation about 3 in the 
power of the velocity dispersion used in the definition of the phase space pseudo- 
density. We suggest that the power 2 may create the more physical quantity, 
whose deviations from a power-law are a diagnostic of incomplete relaxation. 

Subject headings: cosmology: theory — dark matter — halo formation 



1. INTRODUCTION 

Frequently the self-similar in-fall model of structure formation is associated solely with 
the spherically-symmetric, power-law, purely radial, dynamics that was conclusively defined 
in the seminal papers by Fillmore & Goldreich (1984) and Bertschinger (1985). In such 
a restricted formulation, despite the non-linear exactness of the results, this model is not 
considered to have much application to the hierarchical-merging theory of dark matter halo 
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formation. Moreover the model is generally thought to fail to account for the universality of 
form found in the simulated dark matter halos, since it predicts instead a memory of local 
cosmological conditions in the ultimate density power-law. 

However this memory is not in fact very far from universality. If the density of the final 
equilibrium halo, p, is written as a power-law in the radius r in terms of the 'self-similar 
index' (see additional discussion below) j as p oc r~ 2 f , then the index is in fact only weakly 
dependent on 'reasonable' cosmological conditions (Fillmore & Goldreich 1984; Henriksen & 
Widrow 1999) so that there is effective universality. 

This is reflected in the relation for the cosmologically fixed index (Fillmore & Goldreich 
1984; Henriksen & Widrow 1999) j = 3e/(2(e + 1)), where — e is the effective power of the 
initial cosmological peturbation and a 'reasonable' value is close to 2. Our quantity e is 
defined in terms of the density profile while that of Fillmore & Goldreich (1984) is given 
in terms of the mass profile. This means that e = 3ej 9 , where €f g is the parameter used in 
Fillmore & Goldreich (1984). Moreover all strictly radial-orbit simulations show that the 
index j tends universally to 1 if it is initially smaller than 1. This means that the halo 
density profile can never be flatter than r -2 , and it is always close to this value for e « 2. 

As an example of how the self-similar index may be coupled to standard cosmological 
initial conditions, we refer for example to Peacock (1999) (p494) where the linear density 
profile around a maximum in the density field with power spectrum P(k) oc k n gives 2j = 
3(3 + n)/(4 + n) or e = 3 + n. This assumes that the self-similarity is already established in 
this early collapse phase. The estimate works well on the galactic scale where n w — 1, so that 
e 2 as does 2j. In Henriksen & Widrow (1999) the authors suggest instead e = (3 + n)/2, 
which is the expected linear profile around an no peak rather than about a strict maximum 
(Peebles 1993) (p546). This gives better results for cluster scale halos where n ~ 1 and e is 
again nearly 2. It seems that primordial clusters need not dominate their surroundings in 
the same way as do primordial galaxies. 

The strict power-law behaviour is not expected to be maintained at the edge of the 
virialized system (Henriksen & Widrow 1999; Henriksen & Le Delliou 2002; Henriksen 2004). 
This can be due either to mass exhaustion or to tidal truncation (Henriksen 2004, 2006). 
For example in a mass exhaustion situation where e — > oo, the j — > 3/2 and the density 
profile power goes to —3 (Henriksen & Widrow 1999). The tidally truncated edge goes as 
r~ 4 (Henriksen 2004). 

It remains evident despite the preceeding caveats that the cosmologically determined 
self-similarity class can not be extended all the way to the central regions of the halo , 
since the density there has been established as flattened relative to —2j (Navarro, Frenk, 
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& White 1997) ( hereafter NFW), (Taylor & Navarro 2001; Moore et al. 1998). It seems to 
be the strictly radial nature of the orginal self-similar infall models that prevents them from 
experiencing certain kinds of relaxation (MacMillan, Widrow & Henriksen 2006; Huss, Jain, 
& Steinmetz 1999; Barnes et al. 2005a) and hence from evolving a central curved density 
profile. 

It is useful to recall in what follows that the self-similarity index or 'class' is defined more 
generally in Carter & Henriksen (1991) and Henriksen (1997) (see also Henriksen & Widrow 
(1995) for the steady case). In these papers it is shown to give the ratio of the length 
to time dimensions (that is L s /T a ) of whatever quantity is currently controlling the self- 
similarity. The mass dimension is no longer independent because of the need to maintain G 
constant during the scaling(Henriksen & Widrow 1995). This formulation has the advantage 
of allowing the self-similar index, and thus the controlling quantity, to be a parameter in the 
relevant equations; which parameter may change dynamically. 

In Henriksen & Le Delliou (2002) and Henriksen (2004) an asymmetric version of the 
self-similar infall model is treated. In general one can treat anisotropic self-similarity in 
both velocity and real space. Rather than solve the resulting equations by simulation, a 
series expansion is used in terms of the reciprocal resolution in phase space. The smallness 
of this parameter (essentially 1/a with | fixed: see details below) is related to the degree 
of coarse graining both in phase-space and in time. As a — > oo the coarse-graining is 
maximal and the solution is steady. These papers have shown that it is the maximally coarse- 
grained, steady, limit of self-similar infall wherein the density behaves in simple power-law 
fashion. Higher order descriptions of the dynamics reveal a transitory central 'flattening' ( 
always relative to the pure power-law, purely radial profile above) that is presumably due to 
collective relaxation. These terms are not always active however, being notably ineffective 
for the purely radial infall solution. The ultimate central density profile must arise from 
the dynamical evolution of the self-similar index j (since the steady state density power-law 
is just the negative of twice this value, just as for radial orbits discussed above) , which 
phenomenon we refer to as a 'running' index in this paper. 

One might protest that the running index is really a way of saying that the self-similar 
symmetry is lost 'en route' to the relaxed state of the central regions. However simulations 
show an indication of an underlying self-similarity in other ways. One indication is the 
proportional growth of the virial radius and the NFW scale radius (i.e. constant NFW 
'concentration') (MacMillan, Widrow & Henriksen 2006)). Another indication is the power- 
law behaviour of the dimensional or pseudo phase-space density (the real density being the 
distribution function or DF) (Taylor & Navarro 2001; Dehnen & McLaughlin 2005; Barnes 
et al. 2005a, 2006; Austin et al 2005). If then the scaling symmetry is broken during the 
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dynamical relaxation of the central halo, it should only be weakly broken in the sense that 
we can expect a solution in the form of a perturbed self-similarity. If the evolution is slow 
relative to a particle crossing time, then at each epoch the system may be regarded as self- 
similar with the current similarity index. The running index thus corresponds to a kind of 
'adiabatic' self-similarity. It is precisely in this form that we seek a solution for the DF to 
first order. 

So holding the similarity index strictly constant, as is done in the radial self-similar 
infall models, does not allow for effective internal relaxation. There are formal higher order 
terms available for the radial self-similar infall (Henriksen & Le Delliou 2002) which in 
principle describe the approach to a relaxed state, but they do not arise in the absence of 
a suitable mechanism. This has been explicitly confirmed in the simulations by MacMillan, 
Widrow & Henriksen (2006). Analytically, as is emphasized in section 4.1 of this paper, a 
radially biased DF (correct to first order) does not allow a curved density profile to coexist 
with a power-law pseudo- density. Such coexistence is however characteristic of the different 
self-similarity in the halo centre. The purely radial-shell instability detected in Henriksen 
& Widrow (1997) for the original models is thus not the effective mechanism relaxing the 
halo, which is consistent with the moderate redistribution of energy found in these cases 
(Henriksen & Widrow 1999). 

In Henriksen & Widrow (1995) a class of spherically symmetric, stationary, isotropic 
power-law solutions permitting flatter profiles than the cosmologically fixed radial infall 
models was found. These solutions include the Evans & Collett (1997) profile and the NFW 
profile (Henriksen 2006), which the cosmologically fixed radial self-similarity can not do. 
Simulations show (Huss, Jain, & Steinmetz 1999; MacMillan, Widrow & Henriksen 2006) 
that such solutions are closer to the state of the halo centres, although the envelopes are 
closer to the cosmologically fixed radial self-similar infall (Henriksen & Widrow 1999; Lu 
et al. 2006) until edge effects are encountered. Once again a 'running' self-similar index is 
implied, but such steady solutions do not describe the temporal approach to equilibrium 
during which the relaxed state is established. This requires the higher order terms in our 
series expansion. 

As a result of these preliminary considerations, two questions are to be addressed in 
this paper: 

(i) Does the 'similarity class' (the value of the parameter | ) evolve dynamically during 
the halo growth due to collective relaxation in anisotropic systems, and 

(ii) can this evolution be such as to destroy the power-law density profile while main- 
taining a power-law pseudo- density and constant concentration? 
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The answer to the first of these questions was implicitly assumed to be yes in the 
discussion of Henriksen (2006), since there it was argued that the central similarity class of 
an isolated core could be set by conserving particles and phase-space volume rather than by 
the cosmological outer conditions. A 'running' self-similar index is in fact compatible with 
the formulation of self-similarity pioneered in Carter & Henriksen (1991) (see also Henriksen 
(1997) and comment below). 

In this paper therefore we explore the predictions of the first order term in the coarse- 
graining series of the self-similar system (Henriksen & Le Delliou 2002). This term is time 
dependent and includes in a transitory fashion (until it becomes too large) the approach 
to the equilibrium state due to collective relaxation. In addition the dynamical system 
is regarded as being adiabatically self-similar so that the index j is free to vary from a 
value slightly greater than 1 (where it is set by reasonable scale free cosmological conditions 
Fillmore & Goldreich (1984) and above) to a value approaching 0.5 or less as the density 
profile power passes through —1. It is important to understand however that before this 
ultimate steady state is reached, the first order terms can curve the density profile in such a 
fashion as to fit the simulations over a limited radial range. Our conclusion wll be that these 
same terms do not perturb the pseudo phase-space density so strongly in the same radial 
range. 

In the next section after a brief review of the method employed, we calculate the first 
order coarse-grained terms for spherically symmetric self-similarity with velocity anisotropy 
(Henriksen & Widrow 1995; Henriksen & Le Delliou 2002; Henriksen 2004). The spherical 
symmetry should not be thought of as an essential limitation since similar equilibrium solu- 
tions exist with general symmetry (Henriksen 2004) (appendix A). Spherical symmetry does 
however allow additional integrals of the motion which increase the choices available for the 
DF (Kulessa & Lynden-Bell 1992; Henriksen & Widrow 1995; Henriksen 2004). Since there 
are no explicit non-radial forces in spherical symmetry, we are forced to manually repro- 
duce the evolution of the DF from radial orbits to isotropy. The increased choice of DF in 
spherical symmetry permits this. Thus we are able to choose a core DF in section (4.2) that 
corresponds to the spherically averaged equilibrium DF, which form persists even in general 
symmetry (Henriksen 2004). 

Section 3 discusses the calculation of the pseudo phase-space density (referred to fre- 
quently as 'pseudo-density' for brevity) in terms of the higher order DF, and subsequently 
detailed results for several DFs are presented in a section 4. Section five draws our conclu- 
sions and gives a brief speculative discussion. 
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2. Multi-dimensional self-similarity 

2.1. Summary of the Lie symmetry description of self-similarity 

In the mathematical description of self-similarity introduced in Carter & Henriksen 
(1991) and applied to the joint Poisson-Boltzmann system in Henriksen (1997) (with spherical 
symmetry and anisotropic velocity space), one introduces a scaling vector of the form 

b = (a,8,v,X,fi). (1) 

The quantities a, 5, v, A, \i give the logarithmic (additive) form of any scaling in 
the anisotropic phase space of spherically symmetric dynamics, namely, (time, distance, 
velocity, j 2 , mass). The symbol j 2 stands for mean square specific angular momentum. This 
set is complete for dynamical purposes in the sense that the dimension of any physical 
quantity Q may be expressed in terms of its co-vector d^ in this dimension space. For 
example the dimensions of Newton's constant G, namely (in obvious notation) L 3 /(MT 2 ), 
can be expressed as the dimensional co-vector d^ = (—2,3,0,0,-1). This is not a unique 
representation, but since the other dimensions can be expressed in terms of time, length and 
mass; it is the natural one and all other representations are equivalent. 

This description of the self-similar symmetry is that of a Lie symmetry in phase-space 
along the scaling direction k (see Carter & Henriksen (1991) for a more precise mathematical 
description) so that the change in any quantity Q is 

kCQ = (b • d Q )Q, (2) 

and the Lie differentiation is given in spherical symmetry by 

k£ = atd t + 5rd r + vv r d Vr + \j 2 dj2 + fxmd m . (3) 

One proceeds to make the differentiation along k easy by defining a logarithmic time variable 
T satisfying 

d xi T = 0, atd t T = 1, (4) 

and a set of variables, X\ which are invariant under the Lie scaling symmetry so that they 
satisfy 

(k • d x )Xi = 0. (5) 

In these equations x % are the set of physical variables ((r, v r ,j 2 ) in this example) and 
the set X % of the same number are found as characteristic constants of the linear partial 
differential equation (5). The X 1 constitute the self-similar variables and this procedure 
defines 'multi-dimensional' self-similarity. 
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The important feature of the self-similar variables is that when these are used to describe 
the system the scaling operation kCQ becomes simply <9y, and so equation (2) on integration 
yields the multiplicative scaling 

Q = Q (X^) e ^ T . (6) 

The function Q Q is independent of T when the scaling symmetry holds, but if there is no 
such symmetry than it will be in general also a function of T. In such a case the variables 
X 1 are simply a useful set, alternate to the x l . 

In the present example equation (5) gives by the method of characteristics the set X 1 
as (we use notation consistent with previous work) 

R = re~^ aT \ 

Y = v r e~ {u/a)aT , (7) 
Z ee j 2 e-W fl ) aT , 

and the logarithmic time is 

e aT = at. (8) 

There are two additional subtleties of the method. Because for a self-gravitating system 
G is fixed under a scaling operation, we must have b • d G = —2a + 35 — fj, = by equation 
(6). Hence the mass scaling is no longer independent of the space and time scaling, becoming 
rather n = 35 — 2a. The b vector can then be reduced by one dimension. Thus for example 
the mass distribution function f(r,v r ,j 2 ) has dj = (3,-6,0,0,1), which on eliminating \x 
according to the preceeding relation becomes (1,-3,0,0). Then equation (6) allows us to 
write (the factor n is used for convenience) 

nf(r,v r ,j 2 ) ee F(r,v r ,f) = P(R,Y,Z) e -^ aT , (9) 

where P would also depend on T if there were no multi-dimensional self-similarity. 

In like manner, velocity and angular momentum scalings can be written in terms of a 
and 5 as v = 5 — a and A = 45 — 2a. However since there is no invariant velocity in the 
problem v remains free so that 5 ^ a, and similarly 45 ^ 2a since there is no invariant 
angular momentum in the problem and A remains free (near the edge of the system there 
may be tidally induced angular momentum which alters this latter remark (Henriksen 2004), 
as would the presence of a fixed velocity c alter the conclusion for v). 

We see therefore that the parameters a and 5 suffice to express the dimensions of any 
dynamical quantity in our problem, so that from now on we drop the v and A components 
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from b in addition to the [i component as discussed above. Moreover it is only the ratio 
^ or equivalently j = a that appears as a parameter in our transformed equations. This 
parameter a is generally left undetermined so that either appropriate initial conditions, or 
such conditions as arise dynamically, may determine the self-similarity. Following Carter & 
Henriksen (1991) a is referred to as the 'similarity class' since it now determines completely 
the self-similar symmetry. 

The second subtlety is the one introduced in Henriksen & Le Delliou (2002). The phase 
space volume is transformed under the transformation to self-similar variables as 

Ar Av r Af = AR AY AZ e (6/a - 3)aT . (10) 

Thus letting a — * oo gives an ever coarser-graining of the self-similar phase space (a is held 
constant) for fixed increments in the R, Y, Z at fixed T, while in addition giving the temporal 
asymptotic limit by equation (8). Our method of coarse-graining consists in solving the 
transformed equations (expressed in terms of the self-similar variables) in an inverse series 
in a. The zeroth order gives the steady, equilibrium, limit, but the higher order terms 
can describe the approach to this limit in finer detail. By interpreting the coarsest grained 
solution as the equilibrium solution, we are relying on the notion that a fully relaxed system 
should show the same DF at different resolutions, and that this condition is reached after 
an indefinitely long time. 



2.2. Formulation in spherical symmetry with velocity anisotropy 

In this subsection we follow the treatment of Henriksen & Le Delliou (2002) while 
making explicit the units of the various quantities. We recall in our notation that v r is the 
radial velocity and that the specific angular momentum is j 2 = r 2 v\, with v± the transverse 
velocity on a sphere. We also use the logarithmic time of equation (8) and the method of 
the preceeding section. 

It is important to distinguish units from dimensions in this method. We do not assume 
'a priori' essential constants with dimensions of velocity, length, or density, whose existence 
would constrain the scaling symmetry as noted in the previous section. However we do pick 
arbitrary units for these quantities, sometimes called fiducial values. Our various quantities 
are consequently dimensionless numbers with these units understood. 

The unit of the DF is F Q while that of time, radius, velocity and density are r /v , r Q , v Q 
and p respectively. To obtain the equations in the form that we use, these units are related 
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by 

Fo = Po/vl 

(11) 

= 4:iiGp r 2 . 

The potential is measured in units of v 2 . 

The self-similar variables found by the procedure of the previous section are (we repeat 
equation (7) for convenience) 

R = re~^ aT \ Y = v r e~^ aT Z = f e -(^-^T. ( i2 ) 

The remaining dependent variables namely density p and potential $ are expressed according 
to equation (6) using the reduced vector b. These can be written as 

p(r, t) ee Q(R)e- 2aT $(r, t) = T)e 2 ^- l)aT . (13) 

Together with the dependent variable P(R, Y, Z) of equation (9), O and \1/ are the dependent 
variables that express multi-dimensional self-similarity. As remarked in the previous section 
the transformation to these variables does not imply self-similarity unless we set dx = in 
the Boltzmann equation. Otherwise they remain perfectly general, and P, and \& would 
also depend on T. 

The concept of 'running' or adiabatic self-similarity discussed in the introduction is 
compatible with our formulation (that is equation (5) may still be used to find the self- 
similar variables) if S in the equations is replaced by a time dependent quantity 

5 a ee (1/T) J T S(T')dT', (14) 

while the coarse-graining parameter a is left unchanged. Moreover dx acting on P must 
be considered small during a dynamical time. Then one requires 5 a to nearly scale with a 
during the coarse-graining as usual to preserve approximately the similarity class a = a/S a , 
but it must retain sufficient variation in T relative to a to permit j to change by about 
a factor 2 or less between the envelope and core in the application to dark matter halos. 
We shall not in fact need this formulation explicitly in this paper since it suffices simply to 
consider different j in the core and envelope, which procedure corresponds to our similar 
treatment of the evolving DF by fixing the end points. The function S a (T) contains all of 
the internal relaxation physics, and hence is not easily known. 

With these definitions the Poisson equation and the collisionless Boltzmann equation 
(CBE henceforth) become respectively 

±- 2 d R (R 2 d R V) = 0, (15) 
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and 



-d T P (3--l)P + (---R)d R P 

a a a a 



((- - 1)Y + -(d R V - 4)1 d Y P - (4- - 2)Zd z P = 0, (16) 
\ a a Hr J a 



where the scaled density and DF are related by 

1 



= 157 / PdY dZ. (17) 



2.3. Self-similar relaxation 



We proceed as usual by writing a series for P and ^ in inverse powers of a while holding 
the similarity class constant (Henriksen & Le Delliou 2002; Henriksen 2004), which we write 
consistently as a below. We will follow the time dependence of the similarity class only by 
letting a be a variable parameter so that drP = in the equations. Then one finds by 
retaining terms to second order in the governing equations that the solution for P becomes 
(Henriksen & Le Delliou 2002) 

p = P 00 (d, Cl)R a ~' - -Pn(Ci, C 2 2 )iT 3 + 77^22(6, C 2 2 )iT (3+a) (18) 

where the variables Ci and £| are defined as 



Ci = YR^-V, (I = ZR 2(a - 2) (19) 

and are constants on the characteristics of the CBE at all orders of the expansion. The 
functions Pa are only known as solutions to i + 1 equations when the series is terminated at 
the (i + l) th order, as we describe further below. 

The zeroth order (maximally coarse-grained and steady) density and potential become 
respectively (no black hole) using equation (15): 

9 = R~ 2a I 00 , I 00 = J Poo d(i d& (20) 

and (for a ^ 1, which case is logarithmic) 

^ o = - 7oJ R 2 ( 1 -«), 7o = ^ (21) 

'° ' 10 2(a-l)(3-2a) V ; 

Note that the quantities I OQ and 7 Q are dimensionless numbers. 



- 11 - 



As an illustration of the physical implication of these latter relations, we note that 
equations (12), (9) and (13) reveal that to zeroth order nf = P 00 r^ a ~ 3 \ p Q = I 00 r~ 2a and 
$ = — 7 r 2(1 ~ a ). This constitutes a steady solution since Ci and £| are also seen to be 
independent of time, based on the definitions in equation (12). 

The function Pn is given formally by Henriksen & Le Delliou (2002) as corrected in 
Henriksen (2004) as: 

P n = ((a - l)(d 2 - 2 7o ) + C|) d Cl P 00 + 2(a - 2)CiC 2 2 %Poo - Ci(3 - a)P 00 , (22) 
and we find by extending the series to the second order the function P 2 2 to be 

P 22 = ((a - l)(d 2 - 2 7o ) + Qd^Pu + 2(a - 2)GC 2 2 <%Pn - 3CiPn + (3a - 2) 7l <%P 00 . (23) 
The first order potential is given from (15) as (for a ^ 2/3, which case is logarithmic) 

V 1 = _ 7lJ R(2-3a) ^ (24) 

[ 1 3(1 -a)(2-3a)' V 1 

where 

I u EE I P n dd d& (25) 



The nature of the series (18) requires some comment here. It is evident that for any finite 
a the series risks to be non-convergent as R — > 0, depending on the nature of the functions 
Pn and P OQ . Stopping the series at a chosen order by rendering all higher order terms zero 
as in (Henriksen 2004; Henriksen & Le Delliou 2002), will determine these functions (see 
discussion to follow) so as to determine an optimum transitory expression for the DF to the 
required order. However there will always be an inner limit in R to the validity (except as 
a — > oo, which is the completely relaxed steady state) whenever the relevant Pa do not vanish 
with R, and this translates at any r to an upper limit in t. In the non-spherically symmetric 
case, where in principle all of the relaxation physics is included, it may be worth performing a 
renormalization group improvement to the correction term as outlined in McDonald (2006). 

The zeroth order solution is determined in velocity space by requiring P n = (Hen- 
riksen & Le Delliou (2002); Henriksen (2004)), and this yields directly by the method of 
characteristics on what is now equation (22) 

P 00 = K{k)E^), (26) 



where 



o , a—1 



(27) 
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where 

e = \^^- lo \. (28) 

Here K is an arbitrary function of k. 

The special case a = 1 is exponential (Henriksen (2004)) and formally describes the 
singular isothermal sphere, but we do not pursue it here since we believe it to be unstable 
to the development of a core. 

Once again we note that we can pass to physical variables if we wish, by using the 
various definitions and transformations above. These show that 

nf = K(k)\E\ { ^)\ S=\E\r 2{a ~ 1 \ (29) 

and 

\E\ = \£ + £ + K=\E\{f)^\ (30) 

This order is the steady, equilibrium, limit as a — > oo and was originally derived in Henriksen 
& Widrow (1995) although it had been already postulated in a specific form by Kulessa & 
Lynden-Bell (1992). It is interesting that a proper choice of similarity class a allows one to 
include the DF's associated with most commonly used density profiles in their power-law 
limits. Thus referring to the work of Tremaine et al (1994) in the high energy, small r, limit 
(that is for a > 1: it is the low energy, small r limit for a < 1) and identifying their parameter 
T) with our 3 — 2a, one infers from their equation (19) the same DF as in equation (29) for 
3/2 > a > 0. The particular case a = 1/2 yields the central NFW power-law density profile 
and approximate DF. Henriksen (2006) has suggested how this might arise dynamically due 
to conserved phase-space volume in an isolated core. 

In the present work we are concerned with going beyond this limit to the first order 
time-dependent corrections in the DF and the density profile. We will then use the corrected 
DF to calculate the pseudo-density profile in order to test whether this quantity is more or 
less sensitive to the correction. We might hope to find a flattening central density coexisting 
with an accurate power-law in the pseudo phase-space density. 

In order to carry out this plan we need the functions P DO and Pn. By stopping the 
series at first order we require P 22 = 0. Equation (22) and P 22 = are then the expected 
two coupled partial differential equations for the functions P oo (0b Cl) an d Pii(Ci> CD- 

These two equations may be regarded as separate quasi-linear equations with a source 
function that depends on the other unknown function. They have the same characteristics 
in the zeta space namely (£ is measured along a characteristic); 

^ = (a-l)(C 1 2 -2 7o ) + C 2 2 , (31) 
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^ = 2(a-2)C 1 C 2 2 , (32) 

and these yield as a first integral the same characteristic constant as that given in 
equation (27) above. But the characteristic equations for the functions themselves become 

dP 

= Pu + (3 - a)CiP 00 , (33) 

and 

d -w = ^ + W^)^ p ^- (34 > 

The difficulty presents itself in the last term of this latter equation wherein P OQ must be 
known generally in order to evaluate the required derivative on the characteristic in zeta 
space. The only reasonable resolution is by iteration. 

We proceed by substituting the zeroth order expression (26) for P OQ in equation (34) and 
solving for the next approximation to P n by the method of characteristics to find (a ^ 2/3) 

Pu = _ sgn{a _ i) 7li ^)(£)(^)(-J^^ + + K\k){£)^. (35) 

2{a — l) a In k; 

The function K' is arbitrary as is also, we recall, K. One might use this last result to find 
the next order correction for P QO from equation (33), but this will be a small correction and 
will be ignored in what follows. This term could now be transformed into physical variables 
as was done above for the zeroth order, but we shall not need this except to note that unlike 
the zeroth order, the first order is not time independent. By stopping at second order in 
this way the whole series is terminated; but for terms of the same form as the zeroth order, 
which may be absorbed into a renormalized zeroth order term. 

We are now able to calculate directly the density and phase-space pseudo-density be- 
haviour to first order for various choices of the functions K and K', which must be adopted 
on the basis of physical considerations. We proceed to identify the form of these corrections 
in the next sectio and the nature of our choices of DF in the following section. 



3. Density and phase-space pseudo-density to first order 

The density behaviour to first order is readily calculated from equation (17), together 
with the various definitions and the first order series to be 

^ U-^rA. (36) 



I \ cyl 
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This expression shows that the density can 'flatten' relative to its outer slope —2a near 
the centre of the system for two reasons. In the first case (Henriksen & Le Delliou 2002) the 
second term in brackets on the right increases as R — > 0, and over a limited range (see e.g. 
figure 3) the right hand side combines with the power on the left hand side of the equation 
to yield a density profile close to the NFW profile. The active term on the right increases in 
time at fixed r however and so it gives a transitory effect at fixed r, since the density may 
not become negative. Explicitly we must have R a > In/ (al 00 ), which with the definition of 
R becomes r a /t > I n /I 00 . Hence for any fixed a and r there is an upper limit to t for which 
the first order term may be used (see also figure 3 where t is fixed and the minimum r is 
found for validity of the first order). This is true so long as In/I 00 is finite, which is the case 
except for the limiting values a = 1 and a = 2/3. This ratio can be quite small near these 
values where it is essentially linear in a. 

The steady state is found by letting a — > oo, which removes the flattening term on the 
right, and so the second or ultimate flattening must be due to the 'running' nature of the 
self-similarity index a as it decreases through unity to 0.5 or less , as expressed by the power 
on the lhs. The puzzle that we address here is how can the phase-space pseudo-density avoid 
these deviations from a power-law? 

Dchnen & McLaughlin (2005) show that the pseudo-density based on the radial veloc- 
ity dispersion a satisfies as good or better a power-law as that based on the total velocity 
dispersion, so we shall use this for simplicity. In Hansen (2004) Hansen has studied the 
consequences of assuming (both strict and approximate) power-law behaviour for the den- 
sity and the pseudo-density together. He has independently deduced similar results to our 
findings below (for example a = 1/2 is a limiting case) using p/a e for the pseudo- density. 
We define the pseudo-density similarly, but with e = 3e for notational reasons. 

We calculate p/a 3e first rather formally from our various definitions and find (we expect 
e to be of order unity) : 

_P_ r ((2-3e)a+3e) = /gyX 
(T 3e (Cl 2 ) 3e/2 ' 

and so using explicitly our expansion to first order in (17) and the definition of a one obtains 

(l+3e/2) 



3e/2 f 1 - -I^-R- a ) 

P „((2-3e)q+3e) M oo \ aI °° ) /oo\ 

(j3e r(l+3e/2) , „ N 3e/2 " 



( 1 _ M n R-a) 

In this expression the integrals I 00 and I n have already been defined above in eqs. (20) 
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and (25), and we have introduced 

M 00 = J Ci Poo dCi d& (39) 



and 



M n = J Ci Pu d(i d(i (40) 



We observe from eq. (38) that if the right-hand side is roughly constant with R, then 
the power-law predicted for the phase-space pseudo-density is 3e(a — 1) — 2a, that is just 
greater than —2 for a = 1+ and e = 1. As also noted in Hansen (2004), this slope is —2 
independent of a when e = 2/3, and it is approximately —2 independently of e when a ~ 1. 
It is clear that the right-hand side is not strictly constant however, but the question is rather; 
can it be more slowly varying than the correction factor to the density on the right of eq. 
(36)? 

The answer to the previous question requires a tedious calculation of the integrals in 
expression (38) for reasonable choices of the DF, which will be the subject of the next section. 
We note here that if the first order terms in equation (38) are small (they must be less than 
unity), then by first order expansion of the terms in brackets, our calculations must show 
that 

n : r {{2 ^ )a+ ^ - j^r « 1 + ( - ( - 1 )~ A (41) 



o" 3e /£+ 3e / 2 ) al 00 \2 \I U M „ 

is slowly varying compared to the density. This requires the quantity in exterior brackets 
on the right of this latter expression to be small, which we shall refer to as the 'correction 
factor' (C) for brevity, so that 

3e/^Mn \ (42) 

This approximation serves for our present purposes, but there is no reason in principle 
why the more exact expression (38) should not be used. This is so provided that the general 
question is also formulated in such a way that any constants, which depend on particularities 
of the system such as mass, cancel out. This may be achieved by calculating the ratio of the 
logarithmic derivatives of PsD = (p/or 3e )r(( 2 ~ 3e ) a+3e ) and D = pr 2a respectively to obtain 



rlogds = din {PsD) /d In R x din R/ din D = ? ^ 

where we may use equation (42) to write 
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In this approach we would look for a small value of this ratio as a function of a, e and 
r as an idicator of good power law behaviour in the pseudo-density even as D becomes a 
significant function of r. An advantage of this approach is that we do not have the functions 
of I OQ and M OQ on the left hand side of equation (38) to consider, for although these are 
formally constant in r they do depend on a and so implicitly on r as a 'runs'. This is only 
a weak effect when the first order terms are small however and we shall find numerically 
below that ignoring this variation in the linear approach does not yield gross error. It is in 
fact already evident from equation (43) that putting C = renders the exact ratio of the 
logarithmic derivatives as a small quantity of first order. 

4. Distribution Functions and Integrals 

We look at various cases of the DF in order to see when C of equation (41) may indeed 
be small. 



4.1. Envelope-Core transition region 

We know that the orbits in the envelope are biased towards radial orbits while those in 
the core are nearly isotropic (e.g. MacMillan, Widrow & Henriksen (2006)). So we seek to 
describe the transition by choosing a DF that is biased towards radial orbits in the coarse- 
grained limit, but also allows increasing isotropy. We achieve this by choosing in equations 
(26) and (35) respectively 

K = K k(^\ (45) 

and 

K 1 = K 1: (46) 
where K 1 and K Q are constants. This yields 

P 00 = K S^/Cl (47) 

and 

Pn = -^S- 1,2 Kl + K x S^. (48) 
The motivations for this choice of DF are: 

(i) It is evidently consistent with the general coarse-graining series in spherical symmetry 
with velocity anisotropy when stopped at second order, as we have shown above; 
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(ii) The equilibrium term, P OQ , is a radially biased S 1 ^ 2 law Henriksen & Widrow (1999); 

(iii) The first order correction, P n contains a growing isotropic term plus a radially 
biased £~ x l 2 law Fridman & Polyachenko (1984). 

Thus we might expect this DF to describe a transition from primarily radial orbits in a 
steep density profile through the beginning of isotropic flattening. The value of the 'running' 
self-similarity index a should be greater than 1 to describe the steep nature of this transition 
region. 

We shall not give the calculation for this case in detail. It suffices to affirm that we find 
in fact that there is no choice of parameters in this DF which allows the correction factor 
in the pseudo phase-space density to be small while that in the density is substantial. Thus 
we expect both of these quantities to have the power-law behaviour found by setting the 
left-hand sides of equations (41), (36) equal to spatial constants. The self-similar index a 
should then be close to its cosmologically fixed value. 

It is clear that if in fact we do find a region where C can be small while the density is 
substantially flattened (a < 1), then the power of the pseudo- density will steepen in zeroth 
order (see lhs of 41) in this region for e ^ 2/3. We shall see however that C can compensate 
somewhat for this zeroth order behaviour in such a region, as does the limit e — > 2/3. 

We consider in detail the situation in the isotropic, flat, core in the next subsection. 



4.2. Isotropic, Flat, Core 

We turn now to the region that we know to have a flattened density and a power-law 
pseudo-density (Dehnen & McLaughlin 2005; Barnes et al. 2006). We know also that this 
region is close to isotropic (Huss, Jain, & Steinmetz 1999; Barnes et al. 2005a; MacMillan, 
Widrow & Henriksen 2006). Thus we should choose a DF based on eqs.(26) and (35) that 
has a < 1 and is isotropic even in the first order correction term that describes the approach 
to the ultimate steady state. We do this by choosing K = K Q and K' = K\ where both K 
and K\ are constants. Hence in this section we use 

P 00 = K £^h ) (49) 

and 

P n = -| 7l |jf ( — _ a ^ )g(2pi)) + Kl E^>. (50) 
The energy £ is positive in this case, namely £ = (Ci + Cf)/2 + |7o|- 
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We recall our earlier remarks after equation (26) regarding the generality of the zeroth 
order P 00 (S). The term Pu(£) is the first order correction describing the time dependent 
approach to the zeroth order for R not too small. 

It is now straight forward but exceedingly tedious to calculate the quantities I QO , In, 
M OQ and M u . We will discuss only the simplest calculation here to indicate the pattern. But 
accuracy is of the essence and every practical check of the algebra has been made, including 
two independent deductions of the key result. 

The calculation of I QO requires the evaluation of the multiple integral (for brevity we set 
p=(3-a)/(2(a-l))) 



T —IK 

which can be put in the form 



poo poo cp 

/ ^ / d£ h(? 1 n 7 *> (51) 



and hence 



roo poo p 

V2K rfC 2 2 (l7o| + C 2 2 /2)-^ / du^=, (52) 
Jo J i \/u-l 

I 00 = 2V2K ^^B(argl,l/2)\ lo \-^. (53) 
a 



Here B(x,y) is the beta function and 



argl = . (54) 



But by definition (21) (a ^ 2/3) we have I oa = 2(1 — a)(3 — 2a) |^y | so that we obtain a 
relation that is much used in the other results namely 

^ = ^)^w (») 

We note that once the index a is given, the only unknown at zeroth order is K Q . 
Proceeding steadfastly in this fashion we find to this order 

±m±l^\ lo \-^j B (arg2,l/2) 
hi ^ , (ooj 

where 

_ (3-a)(3-2a)aB(arg4,l/2) 
11 3(3a -2)(1- a) B(argl, 1/2)' 1 ' 

We have set 

2 + a 2 - a . . 

ar 9 2 = w —- y argA = —. (58) 
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Moreover provided that a > 1/2 (otherwise the integral over ( 2 must have an arbitrary upper 
cut-off) 

M 00 = 4V2K ^-^\ lo \^B(arg5,3/2), (59) 

and 

argh = . (60) 

Finally we may write using equation (24) 



5a-2 



where 



8^2(1 - a) K x , 
Mu = TT^hol ^B(arg3,3/2) 

OCX 

+ 3 a(3a-2)(l-a) 1701 ^ ^,1,3/2) ( 61 ) 
ar,3 = ^-j. (62) 



All of these preliminaries allow us to write purely as a function of a the ratio that 
appears in C as 

I 00 M U (3-a)(3-2a)(2a-l)B(argl,3/2) 

TdtzO = — — 

IuMoo 3a(3a-2)(l-a) S(ar^5,3/2) 
3(2a - 1) Bjargl, l/2)B(arg3, 3/2) 
5a -2 B(arg2,l/2)B(arg5,3/2) n ' 1 J 

whence follows the explicit correction factor (we repeat equation 44 in a convenient form) 

C = -^-(ratio - 1) - 1. (64) 

The calculation of C depends only on the self-similarity index a and e. We are therefore 
in a position to ask for the function a(e), which is defined by setting C equal to a value 
<< 1. Should this prove to be possible for a range of e about unity, and for values of the 
running index that appear dynamically in the isotropic centre, then we might conclude that 
the pseudo-density with e = 1 has no special significance. Rather it would be the isotropic 
DF, together with the collective relaxation that leads to a < 1 (Huss, Jain, & Steinmetz 
1999; Austin et al 2005; Barnes et al. 2005a; MacMillan, Widrow & Henriksen 2006) that 
produces a power-law pseudo- density for a variety of values of e close to unity (Barnes et al. 
2006). 

We have verified that C is never small according to this type of analysis for an isotropic 
centre with a > 1. Hence such a configuration would require the density and the pseudo- 
density to flatten to the same order relative to their equilibrium power laws, which is not 
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found. This is another indication that a must decrease below unity as the centre of the halo 
is approached 

We shall see in our subsequent discussion that nevertheless (for general e and a running 
a) slight deviations from a pure power-law behaviour for the pseudo- density are predicted 
by our analysis ( there may be a hint of such behaviour in Graham et al. (2005), section 7). 
This is because the steepening relative to —2 produced by the declining index a is not always 
exactly compensated by a flattening due to the first order correction. Barnes et al. (2006) 
however conclude that e = 1 is superior to neighbouring values and that it yields a strict 
power-law greater than —2 . It may be that this result is not yet conclusive (e.g. (Graham 
et al. 2005)), or perhaps as we suggest in the discussion that this is a measure of incomplete 
relaxation in the simulations. This argument requires e = 2/3 to form the ultimate relaxed 
power-law pseudo-density since then the equilibrium power-law of the pseudo-density is 
—2 independent of a ,which can therefore remain fixed near a value that renders C ~ 
throughout the halo. There is as yet no evidence for this state in the simulations. 

In figure (1) C(a, e) is plotted over ranges of a of interest. It is evident that the scenario 
previewed above holds over a rather limited range in a. So we are only permitted to admit 
small variations about the critical a(e) where C = in the figure. 

We see that there is a very precisely defined a such that C = over a range of e close 
to unity. This value is broadened if, for example, we were to be satisfied with C = 0.2. The 
numerical values of 2a required to maintain the pseudo-density power-law lie in the range 
1.4-1.6. 

As a check on our linearized calculation we also present in figure (2) the same calculation 
for the ratio of the logarithmic derivatives at a typical value of x — 0.5. We see that 
numerically the results are very similar, since for e = 1 the ratio of the logarithmic derivatives 
is zero at a = 0.721 compared to C = at a = 0.740. For e = 2/3 the corresponding values 
are 0.773 and 0.800. At x = 1.0 the number for e = 1.0 changes only to 0.728, while at 0.1 
it becomes 0.69. Other values of e yield a similar numerical variation, but do not change 
the results qualitatively. We recall that the major omission from our theory is the ability to 
calculate a(x). This might be found in a detailed fit to a simulation however. 

These density profile powers are associated with the NFW central region (Moore et al. 
1998; Diemand et al. 2005), but they do not correspond to the flattest density profile found 
in the simulations when used in zeroth order (see e.g. eq. 36). However we must remember 
the first order correction factor for the density on the rhs of eq. (36). Indeed the significance 
of our result is that for the appropriate value of a, the density correction factor may be 



Fig. 1. — The figure shows the correction factor C plotted as a function of a and e. The 
horizontal plane corresponds to C — 0. The intersection of the plane with the C(a,e) 
surface defines the curve a(e) for C — 0. The case e = 1 has a zero C for a xs 0.7403, 
while e = 2/3 requires a = 0.8 almost exactly. Intermediate values are e = 0.8, a xs 0.7726; 
e = 1.2,a w 0.7165. 
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0.85 



Fig. 2. — The figure shows the ratio of the logarithmic derivatives given in equation 43 of 
the text namely rlogds(a, e) plotted as a surface over a and e at x = 0.5. The horizontal 
plane shows the zero ratio. The intersection of this plane with the surface defines the curve 
a(e) for zero ratio, which should be compared with the condition C = in figure 1). The 
case e = 1 has a zero ratio for a ~ 0.721, while e = 2/3 requires a = 0.773. Intermediate 
values are e = 0.8. a ps 0.748: e = 1.2. a ps 0.701. 
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substantial while that of the pseudo-density is essentially zero. An example of the fit to the 
NFW density profile that is possible with the form (36) is shown in figure (3) when a = 0.75, 
close to the critical line in figures 1,2). This is only illustrative, but it could be considered a 
reasonable fit to the simulated data of Diemand et al. (2005) (their figure 11) in the range 
10~ 3 r V i r i a i to 10~ 2 r V i r i a i (our x has an arbitrary scale). The figure also shows how sensitive 
the logarithmic slope is to the fitted curve compared to the curve itself. 
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We have noted that the pseudo-density powerlaw should be — ((2 — 3e)a + 3e) —2 in 
terms of the index a in the inner envelope (i.e. near but outside the NFW scale radius) where 
a > 1 Thus there is no reason to suspect that a should vary much from its cosmologically fixed 
value (typically 9/8) until the edge effects near the virial radius destroy the self-similarity. 
We recall that near the virial radius environmental (MacMillan, Widrow & Henriksen (2006)) 
and mass exhaustion effects (Henriksen & Widrow (1999); Lu et al. (2006)) can dominate 
the internal relaxation. As an example of typical outer behaviour we use a = 9/8, which 
gives the pseudo-density power as —(18 — 3e)/8. This gives —15/8 for e = 1 and —2 for 
e = 2/3 (in the latter case the power is —2 strictly independently of a as already remarked). 
These bracket simulated values. 

By way of predictable behaviour inside the scale radius we see from equation (41) that 
the zeroth order pseudo-density steepens (e ^ 2/3) relative to —2 as a decreases from 1, which 
we expect to happen by the arguments in our introduction. This occurs while the density 
is flattening relative to —2 also in zeroth order, although the effect in the pseudo-density is 
weaker as e — > 2/3 . 

The first order correction adds to the steepening for a > a(e) where C > (fig. 1), 
but decreases to zero on the critical curve a(e). After passing through zero, the first-order 
term acts against the zeroth order steepening (figure(l) shows C ~ —0.61 at a = 0.65 for 
e = 1; rlogds 0.57 ), although it is relatively weak relative to the first order density 
correction (the correction factors are not equal in eqs. (36) and (41) or (43) until a = 0.6 
for either e = 1 or e = 0.67). Very similar behaviour is predicted from the exact ratio of the 
logarithmic derivatives as shown in figure (2) provided the density is itself less steep than 

Such a slope oscillation would be detectable in the simulations, but is not strongly 
evident (Dehnen & McLaughlin 2005; Barnes et al. 2006; Graham et al. 2005). The oscillation 
would be less pronounced at smaller e and in fact vanishes at e = 2/3 when the zeroth-order 
power of the pseudo-density becomes independent of a, but then we must remember that C 
is not zero. It vanishes at a = 0.8 for this value of e, but as it is positive at larger a and 
negative at smaller a (figure 1), it would by itself produce an outer steepening followed by 
an inner flattening for this value of e. We must nevertheless emphasize that this remains 
a truer power-law than that of the density, since coincident with this behaviour, equation 
(36) shows that the zeroth and first order terms are combining to flatten the density profile 
monotonically relative to the outer power-law. They can produce the NFW central index of 
— 1 (figure 3) well before a reaches the value 1/2, which would by itself produce this behaviour 
in zeroth order. We do not expect a to be so far from the critical line a(e) however until 
complete relaxation is achieved. 
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Fig. 3. — The left box of the figure shows an example of the fit to the NFW profile (dotted) 
that is possible over a limited region in radius with the first order density formula in the 
text (line). The profiles have been adjusted to have the same slope of —1.1818 and the same 
normalized unit value at x = r/r — 0.1. This requires In/ (al 00 ) ~ 0.053. The fit has been 
made with a = 0.75, which is in the middle of the range that gives zero for the pseudo-density 
correction. One sees that as expected the fit fails at sufficiently small x, but is reasonable 
over a limited range > 0.1. The right box illustrates how sensitive the logarithmic slope is 
to the fitting function over this same range by plotting the logarithmic slope m of the NFW 
function (dotted) and of the fitting function (line). 
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There is one possible exception to the first order slope oscillation expected when e = 2/3 
due to the C variation of figure (l).This will be discussed further in our conclusion section, 
but it requires that a sticks at the value m 0.8 where C vanishes for this e. This would 
produce a rigorous r~ 2 power law for p/a 2 throughout the halo. 

To summarize based on a specific case, the first order flattening of the zeroth order 
power-law a — 3 alluded to above is active when a is in the range ]0.6, 0.74] for e = 1. This 
will tend to preserve the power-law in the pseudo-density which would otherwise steepen in 
zeroth order as a decreases, but there could be detectable steepening as a varies from about 
1 down to 0.74 since C acts also to steepen the pseudo-density in this regime being positive 
(1). However if this transition in the density power happens relatively rapidly in radius, it 
may not be pronounced (e.g. figure 10 (Diemand et al. 2005), appearing only as a small 
'bump' on the curve. Reducing e towards 2/3 will remove the zero-order steepening, but 
makes the first order flattening about equal to that of the density at a = 0.6 ( about one 
half this amount at 0.7, but zero at a = 0.8- all as in figure 1). This should be measurably 
worse than the mean curve for e = 1, if a really does run through this range of values, as is 
concluded in current simulations (Barnes et al. 2006). 

The ratio of the logarithmic derivatives gives us slightly different information than the 
linearized treatment, but it is largely consistent with it. It incorporates the dependence 
on radius that proves to be relatively weak for our current purposes when the first order 
terms are smaller than unity. More importantly figure (2) being roughly independent of 
position, shows that when the density is flatter than r~ 2a on the high a side of the figure, 
the pseudo-density will be steeper than r -(( 2 - 3e ) a + 3e ) anc [ conversely on the low a side. 

5. Discussion and Conclusions 

In this paper we have used two analytic ideas. One was the introduction of the running 
or adiabatic self-similarity, while the other was the series coarse-graining of the self-similarity 
which allows a description of the approach to the final steady state due to internal relaxation. 
It would in principle predict the relaxation in non-spherical symmetry, particularly if a 
renormalization improvement were found. 

We have used the second technique to deduce a DF to zeroth and first order that can 
apply in the flat density, isotropic, limit. Such is expected to be the case inside the NFW 
scaling radius of a halo. This DF appears quite naturally in many places (Tremaine et al 
1994; Henriksen & Widrow 1995) in zeroth order. We used this DF to calculate the first order 
corrections to the density and to the phase-space pseudo-density (generalized slightly by the 
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parameter e). We have found a range of the running index a(e) over which the first order 
correction to the density may be substantial while that of the pseudo-density is negligible. 

More generally, the transitory first order term together with the running index exert 
compensating effects that can maintain the slope of the pseudo-density sensibly constant 
over a wider range of a. This is not expected to change until, and if, a rigorously flat core is 
detected in the simulations. The special case where e = 2/3, so that we consider in fact the 
function p/cr 2 , should show a very good power-law with slope —2 between the virial radius 
and the scale radius, after which the slope will show an oscillation (ultimately flattening ) 
relative to this power-law. 

The region of predominatly radial orbits between the NFW scale radius and the virial 
radius (ignoring edge effects) should be described accurately by the original radial in-fall 
models (Fillmore & Goldreich 1984; Bertschinger 1985), which predict in our terms a slope 
of —((2 — 3e)a + 3e) for the pseudo phase-space density. This is —2 for a — 1, and —15/8 
for a = 9/8 and e = 1. It is also —2 for e = 2/3. 

It is important to note that not just any DF allows the pseudo-density to be a power- 
law while the density flattens. In particular neither a steep density, isotropic DF, nor a 
radially biased, steep density DF does so. Neither the density nor the pseudo-density flattens 
markedly in the envelope where the latter DF is relevant, if edge effects are ignored (e.g. 
(MacMillan, Widrow & Henriksen 2006)). 

Thus we conclude that the remarkable power-law behaviour in the pseudo-density (Tay- 
lor & Navarro 2001) and its related forms is a consequence of the isotropic DF produced in 
the core by processes of relaxation, such as the radial orbit instability (Polyachenko 1981; 
Merritt& Aguilar 1985) . Because the behaviour is found to be similar for a reasonable range 
of e about e = 1 (especially on the low side) we do not think that the form p/a 3 is especially 
significant (see also Hansen (2004), Dehnen & McLaughlin (2005), Barnes et al. (2006)). It 
should be admitted however that currently most simulations (Dehnen & McLaughlin 2005; 
Barnes et al. 2006) do reveal e = 1 as the superior power law although other values of e 
produce power laws superior to that of the density (e.g. (Barnes et al. 2006)) as we would 
expect. Should there prove ultimately to be a preference for e = 1, it would remain a puzzle 
according to the present analysis. It might be attributed to higher order behaviour since 
after all this is a first order calculation. We give an alternate suggestion in our speculation 
below. 

The approximate constancy of the NFW concentration is guaranteed in our treatment 
only by the fact that the evolution is taken to be adiabatically self-similar throughout. 
Together with the previous paragraph, we have thus answered in the affirmative, at least 
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tentatively, the questions posed in the introduction. 

It seems that we have not penetrated very far into the details of the underlying relaxation 
however, beyond suggesting that the equations do know about it! But we note that our 
explanation of the power-law behaviour in the pseudo phase-space density invokes a running 
or adiabatic self-similarity. This is thought to arise as collective effects change the dominant 
conserved quantities from those set by initial conditions to those compatible with an isotropic 
distribution function and an isolated core (Henriksen 2006). The radial orbit instability is 
one collective effect that achieves this (Huss, Jain, & Steinmetz 1999; Barnes et al. 2005a; 
MacMillan, Widrow & Henriksen 2006), but in spherical symmetry we can only imitate such 
effects by using appropriate DF's at the beginning and end of the relaxation process. 

Finally on a speculative note, we observe that pja 2 is predicted here to have a slope 
rigorously —2 until significantly inside the NFW scale radius. Such a relation implies that 
the local Jeans' length (ignoring the difference between the total dispersion and the radial 
dispersion for the moment) scales linearly with radius. Moreover, as discussed in (Henriksen 
2006), we might expect the collective relaxation to proceed by Landau-damping on waves 
created below or at this scale. Hence this relation states that the relaxation is effective up 
to a fixed fraction, probably close to unity, of the local scale. This would be perceived as 
violent 'virialization' on the large scale and asymptotic 'thermalization' on the small scale. 

The flattening deviation from the rigorous inverse square behaviour of this quantity 
would indicate the absence of such relaxation (the Jeans' length becomes greater than the 
local scale). The apparent superiority of e = 1 in the current simulations may be due to the 
expected radial variation of er oc r^ 1_a \ Thus while pja 2 flattens in the absence of complete 
relaxation, dividing by an additional factor a at small r for a < 1 acts against this flattening 
since a is small there. We would predict on these grounds that in simulations with higher 
resolution pjo 2 should become a better power law than p/cr 3 . 

On this view the essential physics is in the behaviour of pjo 2 and the physical running 
value of a would be set by setting C = for e = 2/3. This suggests that the physical 
core value of a would be ~ 0.8 in a perfect simulation. Subsequent flattening of the density 
would be due to first and higher order terms. Fig. (3) is not changed essentially by using 
a = 0.8, so that such flattening is at least possible in a transitory fashion. Ultimately 
the renormalization group approach used by McDonald (2006) may allow a longer lived 
description. 
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